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This paper studies a weakly and asymmetrically coupled three-lane driven diffusive system. A 
non-monotonically changing density profile in the middle lane has been observed. When the extreme value 
of the density profile reaches p=0.5, a bulk induced phase transition occurs which exhibits a shock and a 
continuously and smoothly decreasing density profile which crosses p =0.5 upstream or downstream of the 
shock. The existence of double shocks has also been observed. A mean-field approach has been used to 
interpret the numerical results obtained by Monte Carlo simulations. The current minimization principle 
has excluded the occurrence of two or more bulk induced shocks in the general case of nonzero lane changing 
rates. 
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Driven diffusive system is a rewarding research topic in recent decades, which exhibits non-vanishing 
current even in the steady state and has served as fruitful testing grounds for fundamental research in 
non- equilibrium physics 1 " 5 . The driven diffusive systems exhibit many surprising or counterintuitive 
features, given our experiences with equilibrium systems, e.g., spontaneous symmetry breaking, phase separation, 
etc. 6 " 10 . 

The boundary induced phase transition is another non- equilibrium phenomenon 1114 . For the case of vanishing 
right boundary density, Krug postulated a rather general maximal -current principle that the system tries to 
maximize its stationary current 11 . The maximal -current principle has later been generalized to the extremal 
current principle 1213 



max q(p), p_ >p + 
min q(p), p_ <p + 



(1) 



Here p_ (p+) is the constant effective density of the left (right) reservoir from which particles are flowing into (out 
of) the system and p is system density, which is defined as the occupancy probability of sites. The microscopic 
details of the system only determine the functional form of the current q(p) and the effective boundary densities. 

Motivated by facts such as the unidirectional motion of many motor proteins along cytoskeletal filaments, in 
which motors advance along the filament while attachment and detachment of motors between the cytoplasm 
and the filament occur, the constraint of the conserved dynamics in the bulk in the driven diffusive system has 
been relaxed to consider random particle attachments and detachments in the bulk 15 . The resulting dynamics 
leads to a phase coexistence of low and high density regions separated by a shock. 

The driven diffusive systems have also been extended to include the possibility of transport on multiple parallel 
lanes to describe phenomena such as the extraction of membrane tubes by molecular motors, macroscopic 
clustering phenomena, car traffic and so on (see e.g. Refs. [16-29] and references therein), in which the boundary 
induced phase transitions still could be observed. 

This paper studies a weakly and asymmetrically coupled three-lane totally asymmetric simple exclusion 
process (TASEP), and reports a bulk induced phase transition. Our model is defined in a three-lane lattice of 
NX 3 sites, where N is the length of a lane (see Fig. 1). For each time step, a site is chosen at random. A particle at 
site i and lane j can jump with rate 1 to site i + 1 if it is unoccupied. Otherwise, if site i + 1 is occupied, the particle 
jumps to site i on lane j + 1 with rate oo A and to site i on lane j — 1 with rate oo B (oo B 7^oo A ), provided the target site is 
unoccupied. Obviously, co B =0 for lane 1 and co A =0 for lane 3. We use this rule to mimic the feature of cars, since 
cars usually will not change lane when unhindered. The model was also stimulated by fast and strong motor 
proteins, like kinesins, which do not change lines frequently. However, these details do not change the physics of 
the observed phenomena. Without this rule (i.e., a particle can change lanes provided the target site is empty) the 
system dynamics do not qualitatively change. 
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Figure 1 | Sketch of the three-lane TASEP. The arrow shows allowed 
hopping and the cross shows prohibited hopping. 

The asymmetric coupling makes the middle lane different from 
two other lanes, and bulk induced phase transition occurs on this 
lane. At the boundaries, a particle enters from the left boundary with 
rate oc, and is removed from the right boundary with rate p. As in 
previous studies, a weakly coupling is considered in which co A N=Q A 
and co B N=Q B are kept constant. We have studied the model using 
both mean-field calculation (see section Method) and Monte Carlo 
simulations. 

Results 

Figure 2 shows a typical phase diagram of a three-lane system, in 
which an asymmetrical set of parameters Q A =10 and Q B =0.1 is 
adopted. In the phase diagram, XYZ stands for the state of the three 
lanes. L(H) means that the corresponding lane j is in low (high) 
density phase pj<0.5 (pj>0.5), S means that there is a shock on 
the lane (shocks Si and S 2 denote bulk induced shocks as discussed 
below), D means that there are two shocks on the lane, C means that 
the density profile is a continuously and smoothly decreasing one 
and crosses the density p=0.5. The transitions among LLL, LLS, LSS, 
LLH, LSH, LHH, SSH, SHH, HHH can be easily understood as 
boundary induced phase transitions, see e.g., Refs. [19,20]. For 
instance, when p 3 (l), the density at the right end on lane 3, increases 
to min (P, 1/2) in the LLL phase, a shock is induced from the right 
boundary on this lane and thus the LLL phase transits into LLS phase. 
Fig. 3 shows the typical density profiles of these states. 

We focus on the density profiles in the LLS state. Fig. 4(a) shows 
several typical density profiles with the fixed value P = 0.276. Note 
that in the LLS state, due to existence of a shock on the third lane, the 
density profile is not smooth in the middle lane at the shock location, 
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Figure 2 | Phase diagram of the three-lane TASEP obtained from mean- 
field calculation. The parameters Qa = lOand Q B =0.1. Details about the 
dotted lines and the dashed lines 1-5 are in the text. 



which separates the density profile into two parts. When a is left to 
the dotted line in the phase diagram, the density profile in the down- 
stream part is increasing with x. With the increase of a, the shock 
moves left, so that the downstream part expands. Across the dotted 
line, the density profile in the downstream part becomes non-mono - 
tonically changing. This is because (Ji^ 2 + J3^2) — (/2-»i + J2^3) 
begins to change sign with x in the downstream part. Here J k _^- 
means flow from lane k to lane j. In the LLH state, the shock has 
been expelled out from the left end, and the density profile is always 
monotonically changing. On the boundary between LLH and LSxH, 
the extreme value of the density profile on lane 2 reached p = 0.5. As a 
result, a bulk induced phase transition occurs, which exhibits a shock 
(Si) as well as a continuous and smooth density profile that crosses 
p=0.5 downstream of the shock. Note that when p 2 = 0.5, the coef- 
ficient (l-2p 2 ) in eq.(ll) equals zero. Consequently, the mean field 
equations (10)-(12) cannot be used to solve the density profile in the 
vicinity of p 2 = 0.5, because the numerical solution diverges. We have 
used equations (7)-(9) instead. 

Similarly, in the SHH state, the density profile located in the 
upstream part of the middle lane gradually becomes non-monoton- 
ically changing when (a, p) approaches the corresponding dotted 
line. Across the boundary between LHH and LS 2 H, a bulk induced 
shock (S 2 ) as well as a continuous and smooth density profile that 
crosses p=0.5 upstream of the shock appears, see Fig. 4(b). 

Here we call "bulk induced" phase transition because the shock 
emerges in the bulk. Note that there is another viewpoint, in which 
lanes 1 and 3 are regarded as heterogeneous reservoirs for lane 2. 
Since the heterogeneity on lanes 1 and 3 comes from boundaries, and 
thus the phase transitions are regarded as intrinsically induced by 
boundaries. 

We study the density profiles in the LS X H state. Fig. 5(a) shows 
several density profiles with P = 0.296 is fixed. With the increase of a, 
the shock gradually moves left. Across the boundary between LSiH 
and LCH, the shock is expelled out from the left end, and only the 
continuous and smooth density profile is left in the middle lane. 
Similarly, the shock is expelled out from the right end in the LS 2 H 
phase when across the boundary between LS 2 H and LCH, see 
Fig. 5(b). 

Figs. 6 (a) and (b) show the density profiles in the LSH phase, in 
which (a, p) changes along dashed lines 1 and 2 in the phase diagram, 
respectively. When (a, p) is left of the dotted line, the density profile 
is increasing both upstream and downstream of the shock. However, 
when across the dotted line, the density profile becomes non-mono- 
tonic upstream (downstream) of the shock. On the boundary 
between LSH and LDH, the maximum (minimum) of density profile 
reaches p=0.5. When across the boundary, the bulk induced shock 
appears upstream (downstream) of the first shock and thus two 
shocks exist simultaneously in the middle lane. 

Figs. 6(c) and (d) show the density profiles in the LS X H and LS 2 H 
phase, in which (a, P) changes along dashed lines 3 and 4, respect- 
ively. On approaching the boundary between LS X H (LS 2 H) and LDH, 
P2Q) (p2(0)) gradually approaches P ((1 — a)). Across the boundary, 
a shock is induced from the right (left) boundary, thus, double shocks 
emerge. 

Finally, in the LDH phase, when (a, P) changes along dashed line 
5, the left shock moves toward left and the right shock moves right. 
On the boundary between the LDH phase and the LCH phase, the 
two shocks are expelled from the system, simultaneously, see Fig. 7. 

Discussion 

Now we interpret why only one bulk induced shock is triggered. 
Suppose there are two bulk induced shocks. As a result, the density 
profile crosses p = 0.5 twice, i.e., there are two locations x x and x 2 
(x x <x 2 ) at which p 2 = 0.5. Note that for the double shocks in the LDH 
state, the density profile crosses p=0.5 only once. Thus, the double 
shocks can exist. 
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Figure 3 | Typical density profiles, (a) LLL, (b) LLS, (c) LSS, (d) LLH, (e) LSH, (f) LHH, (g) SSH, (h) SHH, (i) HHH. The scatter data are simulation 
results, and the dashed lines are mean field results, (a, [3) = (0.1, 0.25), (0.1, 0.15), (0.074, 0.087), (0.268, 0.286), (0.2, 0.18), (0.278, 0.166), (0.15, 0.098), 
(0.45, 0.11), (0.45, 0.09) in (a) - (i), respectively. In the simulation, the system size N =20000. 



Now, we have the following equation at the two locations accord- 
ing to the current conservation principle 



- Q A (1 - p 3 ) + 2Q B P 3 2 - Qb(1 - Pi) + 20vPi 2 = 0 



(3) 



Pl(l-Pl) + P 3 ( 1 -P3) = c 



(2) 



where 0 < c < 0.5 stands for a constant. Substituting p 2 = 0.5 into 
eq.(ll), we have 



at locations x x and x 2 . Suppose Q A >Q B , then one has p!<0.5 and 
p 3 >0.5 at the two locations. Thus, from eq.(2), p 3 can be solved 



P3 : 



l + Vl-4[c- Pl (l- Pl )] 



(4) 
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Figure 4 | The density profiles of the middle lane (mean-field results), (a) 

System transits from LLS to LSxH. Parameters are (3 = 0.276, a=0.14, 0.15, 
0.163, 0.17, 0.229, 0.29, 0.308 for curves from bottom to top. (b) System 
transits from SHH to LS 2 H. Parameters are a =0.31, [3=0.088, 0.098, 0.108, 
0.128, 0.152, 0.189, 0.206 for curves from top to bottom. 
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Figure 5 | The density profiles of the middle lane (mean-field results), (a) 

System transits from LS X H to LCH. The parameter |3= 0.296. (b) System 
transits from LS 2 H to LCH. The parameter ot= 0.358. 



SCIENTIFIC REPORTS | 4:5459 | DOI: 1 0.1 038/srep05459 



3 




0.8 



0.0 0.2 0.4 0.6 0.8 1.0 

i/L 



0.0 0.2 0.4 0.6 0.8 1.0 

i/L 
d 

0.81 




0.2 0.4 0.6 0.8 

i/L 



0.2 0.4 0.6 0.8 1.0 

i/L 



Figure 6 | The density profiles of the middle lane (mean-field results), (a) 

shows the transition from LSH to LDH along dashed line 1. Parameters are 
(a, P)= (0.154, 0.173), (0.178, 0.180), (0.190, 0.184), (0.211, 0.190), 
(0.241,0.198) and (0.268, 0.206) for curves from bottom to topatx=0. (b) 
shows the transition along dashed line 2. Parameters are (a, ($)= (0.182, 
0.128), (0.200, 0.145), (0.218, 0.161), (0.234, 0.175) and (0.268, 0.206) for 
curves from top to bottom at x= 1. (c) shows the transition from LSxH to 
LDH along dashed line 3. The parameter a= 0.305. (d) shows the transition 
from LS 2 H to LDH along dashed line 4. The parameter oc=0.318. 



Substituting eq.(4) into (3), one has 
f( Pl )=-0 A < 1 



l + Vl-4[c- Pl (l- Pl ) 



+ 2Q B 



l + Vl-4[c- Pl (l- Pl ) 



(5) 



-Qb(1-Pi) + 2QaP? = 0 

Taking derivative to Pl , we obtain 

df = (l-2 Pl )[Q A + 2Q B (l + d)] 

a Pl d 



+ 4Q aP i+^b 



df 



(6) 



where d= y/l — 4[c — Pl (1 — Pl )] . Since Pl <0.5, j^- >0 is always 

d Pl 

satisfied. This means that eq.(5) has at most one solution. 
Therefore, P i(xi) = P i(x 2 ) and p3 (x 1 ) = P 3(x 2 ). When one lane-chan- 
ging parameter equals 0, both Pl (x) and p3 (x) are monotonically 
changing with x. Therefore, P i(x!) = Pl (x 2 ) and p 3 (x 1 ) = p 3 (x 2 ) can- 
not be satisfied. 

In the general situation that neither Q A nor Q B equals to 0, the 
section x < x x and the section x > x 2 can match each other and 
constitute a shortened system with length 1 — (x 2 — x x ), in which only 
one bulk induced shock exists. The mean-field analysis has shown 
that given the same values of a and p, the flow rate in a shortened 
system is larger than that in the original system, provided there exists 
one bulk induced shock in both systems. Thus, we argue that the 
current minimization principle excludes the occurrence of two or 
more bulk induced shocks. 
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Figure 7 | The density profiles of the middle lane (mean-field results). 

System transits from LDH state to LCH state along dashed line 5. 

Hinsch and Frey have studied a periodic one-dimensional exclu- 
sion process composed of a driven and a diffusive part, and identified 
bulk-driven phase transitions in a mesoscopic limit where both 
dynamics compete 30 . Nevertheless, the system can be regarded as 
two sub -systems connecting together. Therefore, the bulk- driven 
phase transitions are essentially boundary induced ones. 

In our system, we can also treat the location where the density 
profile crosses P = 0.5 as a virtual boundary, which separates the 
system into two sub-systems. For the left sub-system, the effective 
exit rate for the middle lane is P 2 , e ff— 0.5. For the right one, the 
effective entrance rate for the middle lane is a 2 , e ff=0.5. 
Nevertheless, instead of a static boundary, the location of the virtual 
boundary in our system is self-tuned and determined by the values of 
the kinetic rates. 

We also would like to point out that the existence of double shocks 
has been demonstrated when considering detachment and attach- 
ment of particles in the Katz-Lebowitz-Spohn process 31 . However, 
different from the double shocks in LDH state in our model, both 
shocks are boundary induced ones in Ref. [31]. 

Also note that, as pointed out in Ref. [16] , when the minima of the 
current can be realized by several different densities simultaneously, 
multiple shocks can be observed, which separate the density profile 
into several flat segments with different densities corresponding to 
the minima. Nevertheless, the multiple shocks always emerge and 
disappear simultaneously. This is different from the double shocks in 
the present paper. 

This paper only focuses on the three-lane TASEP. When the num- 
ber of lanes further increases, the phase diagram will become much 
more complicated. The questions such as whether two or more bulk 
induced shocks could be observed simultaneously in a lane, whether 
bulk induced shocks can be observed simultaneously in different 
lanes, need to be investigated in the future work. 

Methods 

The hydrodynamic mean field equations of the system could be written as 



at 

dt 



dp 3 

dt ' 



-(1-2 Pi )^-QaPi 2 (1-P 2 ) + ^bP 2 2 (1-Pi) 



-(1 -2p 2 )^ -Q A p 2 2 (l - p 3 ) + Q B p 3 2 (l - p 2 ) 
-QbP 2 2 (1-P 1 ) + ^aP 1 2 (1-P 2 ) 



-(1 -2p 3 )^ -Q B p 3 2 (l - p 2 ) + CW(l - P 3 ) 



(7) 



(8) 



(9) 
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In the steady state, one has 

(1-2 Pi )^- = -QaP 1 2 (1-P 2 )+^bP 2 2 (1-P 1 ) (10) 



(1-2P 2 )-P = -Q A P 2 2 (1-P 3 ) + QbP3 2 (1-P2)-^bP2 2 (1-Pi) + ^AP 1 2 (1-P 2 ) 
dX (11) 



(1 -2p 3 )^ = - P 2 ) + £W(1 - Ps) (12) 

Here pj denotes the density on lane j. Via rescaling the total length to unity (i.e., x=i/ 
L), the boundary conditions pj(0) = min(a, 1/2) and pj(l) = max(l — (3, 1/2) should be 
imposed properly, depending on the specific state of the system. Via numerical 
solving equations (10) - (12) and/or equations (7) - (9), we can obtain the density 
profiles of the three lanes. Specifically, 

In the LLL phase, one has pj(0) = min(a,l/2) (j = 1,2,3). Thus the density profiles of 
the three lanes can be obtained by numerical solving equations (10) - (12) from x=0 
to x=l. 

In the HHH phase, one has pj(l) = max (l-(3, 1/2) (j = 1,2,3). Thus the density 
profiles of the three lanes can be obtained by numerical solving equations (10) - (12) 
from x=l to x=0. 

In the LLS phase, one still has pj(0) = min (a, 1/2). In this case, one needs to assume 
the location of the shock on lane 3, denoted as x s . Then one numerically solves 
equations (10) - (12) from x=0 to x=x s . Denote the density on lane 3 at x=x s as 
p 3 (x s ). One needs to change the density on lane 3 from p 3 (x s ) to 1 — p 3 (x s ) and then 
continues to numerically solve equations (10) - (12) from x=x s to x= 1. If the choice 
of x s is correct, then one would have p 3 (l) = max (1 — (3, 1/2). 

Similarly, in the LSS phase, one still has pj(0) = min (oc, 1/2) and needs to assume the 
locations of the shocks on both lanes 2 and 3 to match p 2 ( 1 ) = p 3 ( 1 ) = max (1 — (3, 1/2). 

In the LLH phase, one has p 1 (0) = p 2 (0) = min (a, 1/2) and needs to assume the 
value of p 3 (0) to match p 3 (l) = max (1 — p, 1/2). 

In the LSH phase, one has p 1 (0) = p 2 (0) = min (a, 1/2) and needs to assume the 
location of the shock on lane 2 and the value of p 3 (0) to match p 2 (l)= p 3 (l)=max 
(1-P, 1/2). 

In the SHH phase, one has pj( 1) = max(l — (3, 1/2) and needs to assume the location 
of the shock on lane 1 to match p 1 (0) = min(a, 1/2). 

In the SSH phase, one still has pj(l)=max(l — (3, 1/2) and needs to assume the 
locations of the shocks on both lanes 1 and 2 to match p 1 (0) = p 2 (0) = min (oc, 1/2). 

IntheLHH phase, one has p 2 (l) = p 3 (l) = max(l — (3, 1/2) and needs to assume the 
value of pi(l) to match p!(0)=min (oc, 1/2). 

In the LCH phase, one has p 1 (0) = min (a, 1/2) and needs to assume the values of 
p 2 (0) and p 3 (0) to match p 3 (l) = max (1 — (3, 1/2). Note that the solution is not unique. 
Due to current minimization principle, the solution corresponding to the minimum 
current is the correct one. 

In the LDH phase, one has pi(0) = p 2 (0) = min (a, 1/2) and needs to assume the two 
locations of the double shocks on lane 2 and the value of p 3 (0) to match p 2 (l) = 
p 3 (l) = max (1 — (3, 1/2). Still the solution is not unique and the current minimization 
principle determines the correct one. 
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